Hola, soy Aitor Gonzalez. Yo aprendo tidymodels para que te sea más fácil de aprender.
Soy data scientist o lo que es lo mismo, estadístico que hace más cosas aparte de una regresión logística. Actualmente trabajo en la unidad de facultad de medicina y biología de la Universidad Autónoma de Barcelona.
Este taller contiene mucho texto. Y suele requerir de un par de veces mirárselo para pillarle todo el tranquillo. No pasa nada si no te enteras a la primera. “Dentro de poco” estará subido a Youtube en mi canal Reestimando.
La historia comienza con este señor de la foto de la ziquierda. Él es Max Khun. Parece un señor muy majo o muy travieso. Yo personalmente creo más lo segundo viendo su creación.
Comienza por aquí por que este señor fue el desarrollador de la librería caret. caret fue una librería que comenzó en 2008 y que asentó las bases del modelado estadístico en R. Desde su creación, el paquete ha sido citado en más de 7000 publicaciones académicas, teniendo un enorme éxito en la comunidad.
El paquete fue un éxito debido a que antes de él, no había una forma de unificar sintaxis entre los principales paquetes de modelado. Cada modelo estaba en un paquete y todos diferentes entre si, un caos para organizar un flujo de trabajo. caret permitía utilizar con una sintaxis común los modelos de varias librerías. Además también traía varias funciones para hacer automáticamente cosas que son necesarias en todo proceso de modelado, como una función que haga una partición de datos en un conjunto de entrenamiento y uno de prueba (training y test).
ELIPSIS en 2012 aparece la el entorno tidyverse de la mano de Hadley Whickam, la otra súper estrella de R en los últimos años. Tidyverse proporciona un montón de herramientas para poder lidiar con facilidad con la manipulación de datos, convirtiéndose en un package súper popular no solo ne R, sino en toda la ciencia datos. El principal punto es la facilidad de sintaxis para todas sus funciones, haciendo que los usuarios tuvieran una lectura de código muy fácil y fluida, enfocándose 100% en la necesidad y no en la programación.
Encantado con un framework en el que toda la sintaxis está ordenada, en 2019, Max se une al tidyverso, y comienza a crear tidymodels cooperando con el equipo de Whickam. El objetivo era convertir Careta algo más modularizado y que permitiera centrar en el modelado en si mismo, más que infinitas líneas de sintaxis para poder crear modelos de machine learning.
A fecha de este taller, noviembre de 2023, hay otra desarrolladora principal de tidymodels además de Max, Julia Silge.
Esta mujer, no solo es un as de la programación, sino que además se dedica a expandir la filosofía tidy a través de sus redes. Su contenido vale oro y diría que es imperdible si quieres no solo aprender más de tidymodels, sino de las posibilidades del aprendizaje automático en general.
- https://www.youtube.com/@JuliaSilge
- https://juliasilge.com/blog/
Entre algunos de sus trabajos más llamativos están uno que me encanta que es Topic modeling for #TidyTuesday Taylor Swift lyrics es decir, “modelizando los temas que tratan las canciones de Taylor Swift” https://juliasilge.com/blog/taylor-swift/
Imaginemos que queremos hacer una regresión lineal. Cogeremos unos datos cualquiera, lo de mtcars.supongamos que queremos hacer una regresión logítica, la sintaxis para eso sería:
Como se pude ver en el ejemplo, la regresión fuera de tidymodels la hacemos con la librería glm. En tidymodels esto se especifica con la función set_engine.
Ahora, veamos que pasaría si queremos hacerla con las liberría h2o
Se va pillando la idea, ¿no?…
tidymodels permite que con cambios mínimos de código puedas cambiar completamente tu flujo de trabajo. Si queremos cambiar un “motor” de modelo, solo tenemos que especificarlo en un solo cambio, y no tener que reescribir todo el código de nuevo. Esto es clave ya que permite jugar mucho a la hora de probar muchos tipos de modelos y ampliar las miras sin aumentar los problemas de código.
Evidentemente, esto no es ni la superficie de tidymodels, aquí solo estamos viendo la isla a lo lejos.
Somos unos expertos en aprendizaje automático y nos piden ayuda de un hospital 🏥
Al parecer, nuestra variable, Var_respuesta, es una enfermedad complicada de diagnosticar. Los médicos no son capaces de determinar si un paciente tendrá o no esa enfermedad al cabo de un año. Entonces han hecho un estudio caso-control en el que se miran unas cuantas variables; y al cabo de un año se da un diagnostico definitivo de var_respuesta. En este caso pues, nos interesa una modelización de tipo clasificatoria binaria, es decir de 2 niveles, o tiene la enfermedad o no la tiene. 👍/👎
Además de de nuestra variable Var_respuesta, tenemos otras 25 variables de diversos tipos, de cognición, bioquímicas, genéticas, etc. El dataframe es una simulación de uno real, e incluye valores perdidos y otras características a tener en cuenta en la fase de pre-procesamiento.
Insisto, el propósito es modelizar un clasifiación binaria
Antes de empezar podemos ver un poco su estructura.
#---------------
# Librerias | Carga de datos ----
#---------------
if (!require('pacman') ) {install.packages('pacman')}
pacman::p_load(
readxl, # Lectura de los datos
tidyverse , # Acceso al entorno de procesamiento "tidyverse"
tidymodels, # Acceso al entorno de modelado "tidymodels"
agua, # Modelaje de entorno h2o
themis, # Soporte de la librería "recipes" para tratar desbalanceo
skimr, # Descriptiva rápida de datos
naniar, # Resumen y Visualización de datos faltantes
knitr, # esta en html para este formato
)
#--- Carga de datos
datos <- read_xlsx('Data/datos.xlsx') %>%
filter(!is.na(Var_respuesta) ) %>%
mutate(Var_respuesta = as_factor(Var_respuesta))
#--- Primeras 10 filas de los datos.
datos %>%
head(.,10) %>%
knitr::kable()| Var_respuesta | Edad | Cognicion_01 | Cognicion_02 | Cognicion_03 | Cognicion_04 | Cognicion_05 | Bioquimica_01 | Bioquimica_02 | Bioquimica_03 | Bioquimica_04 | Bioquimica_05 | Bioquimica_06 | Bioquimica_07 | Bioquimica_08 | Bioquimica_09 | Bioquimica_10 | Bioquimica_11 | Escala_01 | Escala_02 | Escala_03 | Genetica_01 | Genetica_02 | Genetica_03 | Genetica_04 | Genetica_05 | Genetica_06 | Genetica_07 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 14 | 101.82 | NA | NA | NA | NA | NA | NA | 5.65 | 13.65 | 125.71 | NA | NA | 84.29 | 1.48 | 30.35 | 35.68 | 0.00000 | 28.397013 | 7.363183 | NA | 0.9665678 | NA | NA | 1.2921746 | 0.2417377 | NA |
| 0 | 17 | 80.41 | 85.53 | 55.68 | 124.52 | 159.70 | 62.37 | NA | 5.65 | 17.83 | 154.53 | NA | NA | 99.42 | 1.34 | 33.99 | 33.01 | 138.73119 | 16.769985 | 17.764831 | 0.4550571 | NA | NA | NA | NA | NA | 0.0044324 |
| 0 | 14 | 91.00 | 79.76 | 94.28 | 140.68 | 100.07 | 86.33 | 26.24 | 3.18 | 13.09 | 135.88 | 81.53 | 5.75 | 109.45 | 0.95 | 30.44 | 34.21 | NA | 21.382220 | 7.061147 | NA | NA | NA | NA | NA | NA | NA |
| 0 | 21 | 86.99 | 91.33 | 93.65 | 115.00 | 123.01 | 66.84 | 10.23 | 3.48 | 11.97 | 132.15 | 90.66 | 2.58 | 68.60 | 1.04 | 30.10 | 34.00 | 48.83394 | 7.785552 | 20.442332 | NA | NA | NA | NA | NA | NA | NA |
| 1 | 33 | 78.83 | 104.19 | 42.68 | 112.97 | 179.16 | 59.51 | NA | 4.05 | 13.99 | 118.33 | 114.46 | 3.03 | 62.87 | 1.14 | 22.08 | 33.50 | 33.37826 | 28.165315 | 28.817018 | NA | 1.0690024 | 0.0339304 | 0.7820054 | 0.3230982 | NA | 0.3998781 |
| 1 | 27 | 93.54 | 81.16 | 89.32 | 115.83 | 65.07 | 93.17 | 32.35 | 3.50 | 17.95 | 132.61 | 67.61 | 4.10 | 76.30 | 1.11 | 36.50 | 32.24 | 33.45837 | 7.376139 | 22.080631 | 1.3081069 | 0.5763763 | NA | 0.5438819 | NA | 0.4705209 | NA |
| 0 | 16 | 114.25 | NA | NA | NA | NA | NA | 8.14 | 5.68 | 12.25 | 149.89 | NA | NA | 91.22 | 1.28 | 35.86 | 40.11 | NA | 18.489371 | 10.365972 | 1.5581016 | NA | NA | 1.1597700 | NA | 0.6242926 | 0.0579575 |
| 0 | 13 | 74.80 | NA | NA | NA | NA | NA | 93.32 | 4.77 | 8.88 | 164.33 | NA | 4.33 | 84.49 | 1.07 | 23.16 | 27.76 | NA | 18.749210 | 9.068513 | NA | NA | NA | NA | NA | NA | NA |
| 0 | 14 | 66.78 | NA | NA | NA | NA | NA | 54.33 | 2.56 | 16.07 | 161.64 | NA | NA | 99.47 | 0.92 | 30.14 | 34.27 | 133.14928 | 18.116769 | 31.680105 | 1.6247538 | NA | 0.9775707 | NA | 0.9699578 | NA | 1.1983120 |
| 1 | 17 | 84.13 | 100.23 | 69.62 | 137.84 | 91.24 | 78.76 | 28.80 | 4.22 | 14.38 | 151.67 | NA | 5.23 | 91.19 | 1.75 | 18.24 | 37.76 | 113.35402 | 20.327740 | 21.101902 | NA | NA | 1.1801013 | NA | NA | 1.0568568 | NA |
| Name | datos |
| Number of rows | 285 |
| Number of columns | 28 |
| _______________________ | |
| Column type frequency: | |
| factor | 1 |
| numeric | 27 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| Var_respuesta | 0 | 1 | FALSE | 2 | 0: 227, 1: 58 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Edad | 9 | 0.97 | 22.33 | 6.47 | 8.00 | 17.00 | 21.00 | 27.00 | 39.00 | ▂▇▆▅▂ |
| Cognicion_01 | 0 | 1.00 | 86.51 | 29.37 | 0.00 | 78.50 | 90.56 | 103.13 | 143.84 | ▁▁▅▇▂ |
| Cognicion_02 | 82 | 0.71 | 90.22 | 12.98 | 55.23 | 80.74 | 90.43 | 98.98 | 129.15 | ▁▆▇▃▁ |
| Cognicion_03 | 57 | 0.80 | 77.62 | 17.80 | 35.22 | 65.91 | 76.77 | 88.23 | 130.09 | ▂▇▇▃▁ |
| Cognicion_04 | 64 | 0.78 | 134.85 | 51.16 | 15.31 | 97.17 | 137.91 | 172.32 | 245.38 | ▂▅▇▆▂ |
| Cognicion_05 | 70 | 0.75 | 126.23 | 44.38 | 53.89 | 94.13 | 120.96 | 155.32 | 249.33 | ▆▇▆▂▁ |
| Bioquimica_01 | 54 | 0.81 | 75.94 | 14.24 | 42.80 | 66.92 | 74.56 | 84.50 | 118.52 | ▂▇▇▃▁ |
| Bioquimica_02 | 43 | 0.85 | 36.07 | 38.35 | 0.12 | 11.43 | 24.47 | 42.21 | 250.61 | ▇▂▁▁▁ |
| Bioquimica_03 | 13 | 0.95 | 4.05 | 1.64 | 1.27 | 2.90 | 3.67 | 4.74 | 12.68 | ▇▇▂▁▁ |
| Bioquimica_04 | 12 | 0.96 | 14.30 | 1.96 | 8.88 | 12.95 | 14.34 | 15.54 | 20.65 | ▁▆▇▃▁ |
| Bioquimica_05 | 18 | 0.94 | 140.83 | 15.09 | 99.88 | 129.72 | 141.38 | 152.12 | 187.37 | ▁▆▇▃▁ |
| Bioquimica_06 | 54 | 0.81 | 86.32 | 13.75 | 55.96 | 77.21 | 86.04 | 95.10 | 142.74 | ▃▇▅▁▁ |
| Bioquimica_07 | 34 | 0.88 | 3.58 | 1.15 | 0.83 | 3.08 | 3.74 | 4.30 | 6.56 | ▂▂▇▃▁ |
| Bioquimica_08 | 16 | 0.94 | 84.37 | 12.35 | 43.50 | 76.65 | 83.96 | 91.39 | 126.11 | ▁▃▇▃▁ |
| Bioquimica_09 | 57 | 0.80 | 1.13 | 0.29 | 0.51 | 0.91 | 1.12 | 1.33 | 2.09 | ▃▇▆▂▁ |
| Bioquimica_10 | 12 | 0.96 | 30.09 | 3.88 | 18.24 | 27.48 | 30.34 | 32.51 | 40.62 | ▁▃▇▆▁ |
| Bioquimica_11 | 14 | 0.95 | 33.84 | 3.66 | 21.90 | 31.56 | 33.71 | 36.08 | 43.57 | ▁▃▇▆▂ |
| Escala_01 | 32 | 0.89 | 103.18 | 122.93 | 0.00 | 21.23 | 54.69 | 137.77 | 892.17 | ▇▂▁▁▁ |
| Escala_02 | 0 | 1.00 | 18.50 | 8.08 | 6.17 | 11.19 | 18.43 | 23.58 | 41.45 | ▇▇▇▃▁ |
| Escala_03 | 0 | 1.00 | 18.65 | 8.24 | 5.79 | 12.32 | 17.76 | 24.00 | 46.96 | ▇▇▅▂▁ |
| Genetica_01 | 160 | 0.44 | 0.80 | 0.60 | 0.00 | 0.29 | 0.64 | 1.23 | 2.99 | ▇▅▃▁▁ |
| Genetica_02 | 170 | 0.40 | 0.78 | 0.55 | 0.01 | 0.33 | 0.73 | 1.05 | 2.27 | ▇▇▅▂▂ |
| Genetica_03 | 178 | 0.38 | 0.76 | 0.56 | 0.00 | 0.32 | 0.67 | 1.09 | 2.22 | ▇▇▅▂▂ |
| Genetica_04 | 168 | 0.41 | 0.84 | 0.55 | 0.01 | 0.42 | 0.70 | 1.20 | 2.60 | ▇▇▅▂▁ |
| Genetica_05 | 145 | 0.49 | 20.19 | 214.06 | 0.00 | 0.84 | 1.82 | 3.17 | 2534.81 | ▇▁▁▁▁ |
| Genetica_06 | 178 | 0.38 | 0.84 | 0.53 | 0.03 | 0.37 | 0.77 | 1.21 | 2.36 | ▇▇▆▂▁ |
| Genetica_07 | 137 | 0.52 | 0.83 | 0.59 | 0.00 | 0.32 | 0.72 | 1.25 | 2.35 | ▇▆▅▃▂ |
El flujo de trabajo para producir un modelo en tidymodels es el siguiente:
La idea detrás consiste en una sola sintaxis que resuma todo este flujo de trabajo.
Como hemos visto en el ejemplo anterior, con cambiar una sola línea línea, podemos cambiar el tipo de modelo. Y con entender unos pocos puntos clave de todos los nexos, podemos llevar a cabo flujos de trabajo muy eficientes con muy poco coste de mantenimiento y poco cambio de sintaxis.
Para la creación de este “cosmos”, tidymodels utiliza las siguientes librerías:
rsample: Funciones para crear distintos tipos de
remuestreos.recipes: Preprocesador de datos, simplifica el proceso
de preparado de datos para cualquier modelo.parsnip: Interfaz para ajustar modelos, hasta 155 tipos
a fecha de este html.yarstick: Creación de métricas para la evaluación de un
modelo.workflows: Cohesión interna de tidymodels.tune y dials: Creación y gestión de
híper-parámetros de cualquier modelo presente en tidymodels.Para comenzar a crear un modelo de aprendizaje automático es necesario tener datos, fundamental. Estos datos, suelen ser difíciles de conseguir, requieren tiempo, esfuerzo y dinero. Su manejo debe ser optimizado con el fin de abaratar los costes de recopilarles y mejorar los resultados. Además tener en cuenta que el propósito de un modelo de aprendizaje automático en general es que se pueda usar para poder predecir que va a pasar con unos datos que no tenga, o que tenga en un futuro.
Aquí entre el punto del remuestreo. Crearemos un partición de la muestra en unos datos para “entrenar” el modelo y la otra para “evaluarlo”, para ver si tal y como ha quedado calibrado el modelo es útil o no.
Por supuesto, la calidad de esta partición estará ligada a la calidad de los datos. Si los datos recopilados no son suficientemente representativos, el modelo no predecirá bien cuando traigamos casos que no haya visto, aún incluso si la evaluación del modelo ha sido buena.
Básicamente queremos esto:
La librería dentro de tidymodels que nos dará los conjuntos de datos
para hacer esto es rsample.
Algunas de las funciones importantes dentro de este conjunto:
initial_split(): Principal función dentro de la
librería.Crea un objeto “split” sencillo para hacer una
partición de datos# 🔴 La función principal de Rsample, hace una partición de un conjunto de 🔴
initial_split(datos, prop = 0.6 )## <Training/Testing/Total>
## <171/114/285>
bootstraps(): crea un tibble con n cantidad de objetos
“split” y su correspondiente identificador. Estas muestras son
de tipo bootstrap, es decir, muestras en las que un individuos puede
aparecer más de una vez en el conjunto de datos, tanto en entrenamiento
como en evaluación.## # Bootstrap sampling
## # A tibble: 5 × 2
## splits id
## <list> <chr>
## 1 <split [285/99]> Bootstrap1
## 2 <split [285/103]> Bootstrap2
## 3 <split [285/113]> Bootstrap3
## 4 <split [285/104]> Bootstrap4
## 5 <split [285/113]> Bootstrap5
rolling_origin(): Crea un tibble con “n”
cantidad de objetos “split” y su correspondiente identificador.
En este caso, los remuestreos no son aleatorios y cada partición
contienen datos que consecutivos. Es decir, se comienza por la fila 1 y
se hacen bloqueen a cantidades iguales. Se usa sobretodo en la partición
de datos de tiempo o cuando se está estudiando la creación de una
partición optimizada entre el conjunto de entrenamiento y de
evaluación.## # Rolling origin forecast resampling
## # A tibble: 6 × 2
## splits id
## <list> <chr>
## 1 <split [20/20]> Slice1
## 2 <split [20/20]> Slice2
## 3 <split [20/20]> Slice3
## 4 <split [20/20]> Slice4
## 5 <split [20/20]> Slice5
## 6 <split [20/20]> Slice6
ggpubr::ggarrange(
rolling_origin(datos, initial = 20, assess = 20, skip = 40, cumulative = FALSE) %>%
tidy() %>%
ggplot(aes(x = Resample, y = factor(Row), fill = Data)) +
geom_tile() +
scale_fill_manual( values = c('cyan','orange'))
,
rolling_origin(datos, initial = 50, assess = 20, skip = 70, cumulative = FALSE) %>%
tidy() %>%
ggplot(aes(x = Resample, y = factor(Row), fill = Data)) +
geom_tile() +
scale_fill_manual( values = c('cyan','orange'))
,
rolling_origin(datos, initial = 20, assess = 20, skip = 40, cumulative = T) %>%
tidy() %>%
ggplot(aes(x = Resample, y = factor(Row), fill = Data)) +
geom_tile() +
scale_fill_manual( values = c('cyan','orange'))
,
rolling_origin(datos, initial = 50, assess = 50, skip = 70, cumulative = T) %>%
tidy() %>%
ggplot(aes(x = Resample, y = factor(Row), fill = Data)) +
geom_tile() +
scale_fill_manual( values = c('cyan','orange'))
,
common.legend = T, legend = 'bottom'
)Una de las cosas que pasan habitual a la hora de crear una partición de datos es que, al ser aleatorias, uno de los conjuntos contenga más datos de un tiempo de otro tipo.
Idealmente querríamos que el conjunto de entrenamiento y que el conjunto de evaluación contuvieran la misma proporción de las categorías.
Para ello, en varias de las funciones de rsample se cuenta con el argumento strata que permite hacer esto de forma estratificada para una variable o variables de interés, es este caso, para nuestros datos y nuestra variable conjunta:
set.seed(4147)
# Splits ----
## Sin estratificar
Split_datos_no_strat <- initial_split(datos, prop = 0.70 )
# Put 3/4 of the data into the training set
# 🔴 Initial_split con datos estratificados 🔴
## Estratificados
Split_datos_strat <- initial_split(datos, prop = 0.70, strata = Var_respuesta)ggpubr::ggarrange(
ggpubr::ggarrange(
training(Split_datos_no_strat) %>%
group_by(Var_respuesta) %>%
summarise(recuento = n()) %>%
ungroup() %>%
mutate(prop_var_respuesta = recuento / sum(recuento) ) %>%
ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
geom_bar(stat = "identity") +
geom_text(aes(label = scales::percent(prop_var_respuesta)), position = position_stack(vjust = 0.5), size = 4) +
scale_fill_manual( values = c("azure4", "azure3"))
,
testing(Split_datos_no_strat) %>%
group_by(Var_respuesta) %>%
summarise(recuento = n()) %>%
mutate(prop_var_respuesta = recuento / sum(recuento) ) %>%
ungroup() %>%
ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
geom_bar(stat = "identity") +
geom_text(aes(label = scales::percent(prop_var_respuesta)), position = position_stack(vjust = 0.5), size = 4) +
scale_fill_manual( values = c("azure4", "azure3")),
common.legend = T, legend = 'bottom'
),
ggpubr::ggarrange(
training(Split_datos_strat) %>%
group_by(Var_respuesta) %>%
summarise(recuento = n()) %>%
ungroup() %>%
mutate(prop_var_respuesta = recuento / sum(recuento) ) %>%
ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
geom_bar(stat = "identity") +
geom_text(aes(label = scales::percent(prop_var_respuesta)), position = position_stack(vjust = 0.5), size = 4) +
scale_fill_manual( values = c("#E86EB7", "#108064"))
,
testing(Split_datos_strat) %>%
group_by(Var_respuesta) %>%
summarise(recuento = n()) %>%
mutate(prop_var_respuesta = recuento / sum(recuento) ) %>%
ungroup() %>%
ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
geom_bar(stat = "identity") +
geom_text(aes(label = scales::percent(prop_var_respuesta)), position = position_stack(vjust = 0.5), size = 4) +
scale_fill_manual( values = c("#E86EB7", "#108064")),
common.legend = T,
legend = 'bottom'
),
labels = c('No estratificados', 'Estratificados'), label.x = 0.1
) %>%
plot()| Aclaración: los colores elegidos para el gráfico son de mi escena preferida de Spiderman, into the spiderverse. dejad siempre “huevos de Pascua” en los trabajos, hace que las cosas sean más emocionantes :P |
Para esto usaremos las funciones training() y test()
# 🔴 Deshacemos las particiones con las funciones training y testing 🔴
Training_datos <- training(Split_datos_strat)
Testing_datos <- testing(Split_datos_strat)
knitr::kable(head(Training_datos,10))| Var_respuesta | Edad | Cognicion_01 | Cognicion_02 | Cognicion_03 | Cognicion_04 | Cognicion_05 | Bioquimica_01 | Bioquimica_02 | Bioquimica_03 | Bioquimica_04 | Bioquimica_05 | Bioquimica_06 | Bioquimica_07 | Bioquimica_08 | Bioquimica_09 | Bioquimica_10 | Bioquimica_11 | Escala_01 | Escala_02 | Escala_03 | Genetica_01 | Genetica_02 | Genetica_03 | Genetica_04 | Genetica_05 | Genetica_06 | Genetica_07 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 14 | 66.78 | NA | NA | NA | NA | NA | 54.33 | 2.56 | 16.07 | 161.64 | NA | NA | 99.47 | 0.92 | 30.14 | 34.27 | 133.149284 | 18.116769 | 31.680105 | 1.6247538 | NA | 0.9775707 | NA | 0.9699578 | NA | 1.1983120 |
| 0 | 23 | 108.18 | 84.49 | 82.39 | 183.52 | 166.31 | 101.24 | 3.66 | 3.18 | 20.29 | 148.58 | 100.23 | 3.53 | 81.72 | 1.11 | 35.75 | 34.03 | 68.943373 | 12.183308 | 15.670177 | NA | 1.797438 | NA | NA | 0.1373531 | 1.1481697 | NA |
| 0 | 26 | 0.00 | NA | NA | NA | NA | NA | 129.41 | 5.99 | 14.98 | 145.04 | 94.53 | 4.20 | 78.76 | 0.87 | 32.57 | 26.89 | 2.898493 | 11.236520 | 6.628174 | NA | NA | NA | 1.6622978 | 2.6878823 | 1.6836672 | 0.1001570 |
| 0 | 14 | 52.88 | NA | NA | NA | NA | NA | NA | 2.06 | 13.98 | 132.94 | 105.75 | NA | 82.49 | NA | 34.73 | 29.83 | 99.600374 | 25.026820 | 12.949102 | 0.1247104 | NA | NA | 0.6977428 | 1.2486785 | 0.4195379 | 0.3760827 |
| 0 | 27 | 109.25 | 88.04 | 103.31 | 143.93 | 144.97 | 100.08 | 8.09 | 4.49 | 19.09 | 187.37 | 84.47 | 3.87 | 91.97 | 1.26 | 30.74 | 34.95 | 0.000000 | 16.798451 | 7.607729 | 0.8493442 | NA | 0.8781136 | 1.9962258 | NA | NA | NA |
| 0 | 30 | 105.60 | 104.09 | 86.45 | 149.99 | 105.02 | 93.65 | 7.20 | 4.51 | 15.65 | 145.00 | 90.13 | 3.57 | 90.67 | 1.05 | 28.98 | 31.50 | 51.087736 | 6.731747 | 10.894387 | NA | 1.349131 | 0.3128530 | NA | 1.9916249 | 1.2268437 | NA |
| 0 | 12 | 91.73 | NA | NA | NA | NA | NA | 30.29 | 3.15 | 13.94 | 134.92 | 103.40 | 3.97 | 96.92 | 1.27 | 21.42 | 29.26 | 160.752695 | 17.851526 | 18.519326 | 0.8739840 | NA | 0.5427362 | NA | 1.2463796 | 0.8887141 | 0.7917460 |
| 0 | 17 | 91.97 | 95.35 | 80.65 | 172.32 | 122.78 | 68.50 | 40.55 | 2.77 | 12.49 | 145.84 | 74.00 | 4.01 | 86.12 | 1.16 | 27.01 | 33.90 | 15.319511 | 14.058768 | 17.484487 | 0.2352628 | NA | 0.3039570 | NA | NA | NA | 0.9261596 |
| 0 | 14 | 69.65 | NA | NA | NA | NA | NA | 55.17 | 4.87 | 11.98 | 165.35 | 91.39 | 5.73 | 52.67 | 1.21 | 29.25 | 32.46 | 76.930164 | 22.999587 | 14.681546 | NA | NA | NA | NA | NA | NA | NA |
| 0 | 34 | 111.50 | 58.99 | 101.17 | 181.40 | 119.87 | 106.90 | 9.96 | 3.19 | 14.47 | 159.95 | 102.97 | 3.98 | 85.87 | 0.91 | 30.43 | 29.18 | 7.373436 | 8.195322 | 6.718927 | 1.2862183 | NA | NA | NA | NA | NA | 0.2123626 |
| Var_respuesta | Edad | Cognicion_01 | Cognicion_02 | Cognicion_03 | Cognicion_04 | Cognicion_05 | Bioquimica_01 | Bioquimica_02 | Bioquimica_03 | Bioquimica_04 | Bioquimica_05 | Bioquimica_06 | Bioquimica_07 | Bioquimica_08 | Bioquimica_09 | Bioquimica_10 | Bioquimica_11 | Escala_01 | Escala_02 | Escala_03 | Genetica_01 | Genetica_02 | Genetica_03 | Genetica_04 | Genetica_05 | Genetica_06 | Genetica_07 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 14 | 101.82 | NA | NA | NA | NA | NA | NA | 5.65 | 13.65 | 125.71 | NA | NA | 84.29 | 1.48 | 30.35 | 35.68 | 0.00000 | 28.397013 | 7.363183 | NA | 0.9665678 | NA | NA | 1.292175 | 0.2417377 | NA |
| 0 | 17 | 80.41 | 85.53 | 55.68 | 124.52 | 159.70 | 62.37 | NA | 5.65 | 17.83 | 154.53 | NA | NA | 99.42 | 1.34 | 33.99 | 33.01 | 138.73119 | 16.769985 | 17.764831 | 0.4550571 | NA | NA | NA | NA | NA | 0.0044324 |
| 0 | 14 | 91.00 | 79.76 | 94.28 | 140.68 | 100.07 | 86.33 | 26.24 | 3.18 | 13.09 | 135.88 | 81.53 | 5.75 | 109.45 | 0.95 | 30.44 | 34.21 | NA | 21.382220 | 7.061147 | NA | NA | NA | NA | NA | NA | NA |
| 0 | 21 | 86.99 | 91.33 | 93.65 | 115.00 | 123.01 | 66.84 | 10.23 | 3.48 | 11.97 | 132.15 | 90.66 | 2.58 | 68.60 | 1.04 | 30.10 | 34.00 | 48.83394 | 7.785552 | 20.442332 | NA | NA | NA | NA | NA | NA | NA |
| 0 | 16 | 114.25 | NA | NA | NA | NA | NA | 8.14 | 5.68 | 12.25 | 149.89 | NA | NA | 91.22 | 1.28 | 35.86 | 40.11 | NA | 18.489371 | 10.365972 | 1.5581016 | NA | NA | 1.159770 | NA | 0.6242926 | 0.0579575 |
| 0 | 13 | 74.80 | NA | NA | NA | NA | NA | 93.32 | 4.77 | 8.88 | 164.33 | NA | 4.33 | 84.49 | 1.07 | 23.16 | 27.76 | NA | 18.749210 | 9.068513 | NA | NA | NA | NA | NA | NA | NA |
| 1 | 17 | 84.13 | 100.23 | 69.62 | 137.84 | 91.24 | 78.76 | 28.80 | 4.22 | 14.38 | 151.67 | NA | 5.23 | 91.19 | 1.75 | 18.24 | 37.76 | 113.35402 | 20.327740 | 21.101902 | NA | NA | 1.180101 | NA | NA | 1.0568568 | NA |
| 0 | 31 | 128.68 | NA | 99.23 | 173.48 | 144.55 | 101.62 | 21.34 | 5.61 | 14.16 | 147.96 | 104.95 | 4.06 | 91.39 | 1.23 | 31.89 | 32.58 | 118.10758 | 17.200513 | 14.212995 | NA | NA | NA | 0.972422 | 1.642294 | NA | 0.5527762 |
| 0 | 20 | 110.69 | 88.69 | 72.28 | 143.15 | 193.88 | 86.83 | 25.79 | 4.50 | 13.12 | 134.66 | 68.58 | 3.56 | 81.48 | 1.11 | 28.27 | 31.64 | 225.00452 | 17.789336 | 20.238253 | NA | NA | NA | NA | NA | NA | NA |
| 0 | 16 | 83.73 | NA | NA | NA | NA | NA | 23.32 | 4.16 | 12.12 | 137.10 | 83.00 | 5.53 | 71.56 | 0.99 | 32.46 | 35.09 | 133.40626 | 20.677244 | 15.434748 | NA | NA | NA | NA | NA | NA | NA |
parsnipEl pre-procesamiento es clave en el proceso de creación de un modelos de machine learning. Cuando queremos hacer un dataframe para modelizar en aprendizaje automático, no viene todo arreglado. Incluso con una recopilación de datos buena puede que necesitemos algunos arreglos. Algunas de las causas de esto pueden ser:
Hacer todo estos procesos puede ser engorroso, pero podemos salir del paso fácilmente gracias por ejemplo de dplyr. Sin embargo ¿Cómo hacemos cuando queramos poner un modelo a grane scala y entren datos a cada momento? ¿Debemos ejecutar el pre-procesamiento una y otra vez?.
La clave para resolver esto está en el paquete recipes. Esta librería trae consigo todo un conjunto de funciones que permiten un pre-procesamiento fino, fino.
Para poder elaborar un flujo completo de pre-procesamiento debemos seguir los siguuientes pasos:
recipe(): la base. Requiere de qué fórmula queremos
para el modelo y de qué conjunto de datos partimos. Los datos aquí
incorporados serán los datos de entrenamientostep(): el proceso que queramos aplicar a nuestros
datos, cualquier transformación de datos va aquí.prep(): la receta se prepara, es decir se aplica al
conjunto de entrenamiento y se abstrae para un modelado cualquierabake(): pone el pre-procesamiento en producción para un
conjunto de datos cualquiera.Hagamos un ejemplo básico: queremos una receta que sencillamente impute los datos faltantes por la media correspondiente a cada variable. Para eso vamos con ejemplo sencillo. La siguiente receta dejará los datos preparados para un modelo de aprendizaje automático, pero sin aplicarle ninguna transformación.
Training_datos %>%
# 🔴 1. Recipe (crea la formula generalizada del modelo)
recipe(
as.formula("Var_respuesta ~ Cognicion_01 +Cognicion_02+Escala_01+Escala_02"),
data = . ,
strata = Var_respuesta ) %>%
# 🔴 2. step
# 🟡 En este caso no aplicamos ninguina trasnformación
# Step_log
# 🔴 3. Prep (asienta la receta y generalizala)
prep() %>%
# 🔴 4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
bake(., new_data = NULL) %>%
# mostramos resultados
head(.,10) %>%
knitr::kable()| Cognicion_01 | Cognicion_02 | Escala_01 | Escala_02 | Var_respuesta |
|---|---|---|---|---|
| 66.78 | NA | 133.149284 | 18.116769 | 0 |
| 108.18 | 84.49 | 68.943373 | 12.183308 | 0 |
| 0.00 | NA | 2.898493 | 11.236520 | 0 |
| 52.88 | NA | 99.600374 | 25.026820 | 0 |
| 109.25 | 88.04 | 0.000000 | 16.798451 | 0 |
| 105.60 | 104.09 | 51.087736 | 6.731747 | 0 |
| 91.73 | NA | 160.752695 | 17.851526 | 0 |
| 91.97 | 95.35 | 15.319511 | 14.058768 | 0 |
| 69.65 | NA | 76.930164 | 22.999587 | 0 |
| 111.50 | 58.99 | 7.373436 | 8.195322 | 0 |
Training_datos %>%
# 🔴 1. Recipe (crea la formula generalizada del modelo)
recipe(
as.formula("Var_respuesta ~ Cognicion_01 +Cognicion_02+Escala_01+Escala_02"),
data = . ,
strata = Var_respuesta ) %>%
# 🔴 2. Step (haz la transofrmación que requieras)
# 🟡 Este step será el logaritmo de todas las variables predictoras númericas.
step_log(all_numeric(), -all_outcomes()) %>%
# 🔴 3. Prep (asienta la receta y generalizala)
prep() %>%
# 🔴 4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
bake(., new_data = NULL) %>%
# mostramos resultados
head(.,10) %>%
knitr::kable()| Cognicion_01 | Cognicion_02 | Escala_01 | Escala_02 | Var_respuesta |
|---|---|---|---|---|
| 4.201404 | NA | 4.891471 | 2.896838 | 0 |
| 4.683796 | 4.436633 | 4.233286 | 2.500067 | 0 |
| -Inf | NA | 1.064191 | 2.419169 | 0 |
| 3.968025 | NA | 4.601166 | 3.219948 | 0 |
| 4.693639 | 4.477791 | -Inf | 2.821287 | 0 |
| 4.659658 | 4.645256 | 3.933544 | 1.906835 | 0 |
| 4.518849 | NA | 5.079867 | 2.882089 | 0 |
| 4.521462 | 4.557554 | 2.729127 | 2.643246 | 0 |
| 4.243483 | NA | 4.342898 | 3.135476 | 0 |
| 4.714025 | 4.077368 | 1.997884 | 2.103563 | 0 |
¡WALA! Datos pre-procesados y listos.
A partir de aquí en cielo es el límite ya que hay steps para todo, recomiendo leer la librería recipes y ver que tienen disponible. Al final de esta sección dejaremos lista una receta como ejemplo
Hacer una imputación simple suele no ser la mejor vía de hacerlo. lo ideal siempre es imputación múltiple. Es decir, una imputación que rellena los valores que faltan con números plausibles derivados de las distribuciones y relaciones entre las variables observadas en el conjunto de datos.
En tidymodels se puede hacer imputación múltiple, y valga la redundancia, de múltiples formas. Por ejemplo, podemos imputar una variable usando una regresión lineal a partir de otras. Hagamos una prueba con la variable Cognicion_02 a partir de Cognicion_01 y Escala_02
Training_datos %>%
# 🔴 1. Recipe (crea la formula generalizada del modelo)
recipe(
as.formula("Var_respuesta ~ Cognicion_01 +Cognicion_02+Escala_01+Escala_02"),
data = . ,
strata = Var_respuesta ) %>%
# 🔴 2. Step (haz la transofrmación que requieras)
# 🟡 Este step será la imputación por regresiópn lineal con otras covariables
step_impute_linear( Cognicion_02, impute_with = imp_vars(Cognicion_01, Escala_02) ) %>%
# 🔴 3. Prep (asienta la receta y generalizala)
prep() %>%
# 🔴 4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
bake(., new_data = NULL) %>%
# mostramos resultados
head(.,10) %>%
knitr::kable()| Cognicion_01 | Cognicion_02 | Escala_01 | Escala_02 | Var_respuesta |
|---|---|---|---|---|
| 66.78 | 95.30378 | 133.149284 | 18.116769 | 0 |
| 108.18 | 84.49000 | 68.943373 | 12.183308 | 0 |
| 0.00 | 105.71896 | 2.898493 | 11.236520 | 0 |
| 52.88 | 97.79188 | 99.600374 | 25.026820 | 0 |
| 109.25 | 88.04000 | 0.000000 | 16.798451 | 0 |
| 105.60 | 104.09000 | 51.087736 | 6.731747 | 0 |
| 91.73 | 91.30365 | 160.752695 | 17.851526 | 0 |
| 91.97 | 95.35000 | 15.319511 | 14.058768 | 0 |
| 69.65 | 95.03224 | 76.930164 | 22.999587 | 0 |
| 111.50 | 58.99000 | 7.373436 | 8.195322 | 0 |
Esto va genial cuando conocemos la estructura de los datos y sabemos con certeza la dependencia de unas variables con otras para imputarlas fácilmente, en lugar de confiar ciegamente en algoritmos de imputación.
Como parte del pre-procesamiento, a veces puede ser necesario convertir ciertas variables a PCA para usarlas como componentes. No entraré en su uso, solo especificar la forma de hacerlo.
Training_datos %>%
# 🔴 1. Recipe (crea la formula generalizada del modelo)
recipe(
as.formula("Var_respuesta ~ Cognicion_01 + Cognicion_02 + Cognicion_03 + Escala_01 + Escala_02 + Escala_03"),
data = . ,
strata = Var_respuesta ) %>%
# 🔴 2. Step (haz la transofrmación que requieras)
# 🟡 Este step será la imputación por regresiópn lineal mcon otras covariables
step_impute_bag(all_numeric()) %>%
step_pca( Cognicion_01 , Cognicion_02 , Cognicion_03, num_comp = 2, id = "pca") %>%
# 🔴 3. Prep (asienta la receta y generalizala)
prep() %>%
# 🔴 4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
tidy(id = "pca") %>%
# mostramos resultados
head(.,10) %>%
knitr::kable()| terms | value | component | id |
|---|---|---|---|
| Cognicion_01 | -0.5964794 | PC1 | pca |
| Cognicion_02 | -0.6133048 | PC1 | pca |
| Cognicion_03 | -0.5177543 | PC1 | pca |
| Cognicion_01 | 0.7421220 | PC2 | pca |
| Cognicion_02 | -0.6671338 | PC2 | pca |
| Cognicion_03 | -0.0647105 | PC2 | pca |
| Cognicion_01 | 0.3057241 | PC3 | pca |
| Cognicion_02 | 0.4228354 | PC3 | pca |
| Cognicion_03 | -0.8530785 | PC3 | pca |
ADASYN, o Adaptive Synthetic Sampling, es un algoritmo bastante utilizado cuando tenemos conjuntos de datos desequilibrados.
Los conjuntos de datos desequilibrados se producen cuando la distribución de clases en los datos es desigual, lo que dificulta el entrenamiento de modelos que pueden estar sesgados hacia la clase mayoritaria.
# 🔴 Libreria que permite instalar algoritmos e sobremuestreo en tidymodels.
pacman::p_load(themis)
Training_datos_adasyn <- Training_datos %>%
# 1. Recipe (crea la formula generalizada del modelo)
recipe(
as.formula("Var_respuesta ~ Cognicion_01 +Cognicion_02+Escala_01+Escala_02"),
data = . ,
strata = Var_respuesta ) %>%
# 2. Step (haz la transofrmación que requieras)
# step de imputación por arboles manteniendo la estrucura de todas las covariables
step_impute_bag(all_numeric()) %>%
# 🔴 step para sobremuestrear la el nivel menor en la variable respuesta.
step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>%
# 3. Prep (asienta la receta y generalizala)
prep() %>%
# 4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
bake(., new_data = NULL)| Cognicion_01 | Cognicion_02 | Escala_01 | Escala_02 | Var_respuesta |
|---|---|---|---|---|
| 66.78 | 102.05505 | 133.149284 | 18.116769 | 0 |
| 108.18 | 84.49000 | 68.943373 | 12.183308 | 0 |
| 0.00 | 99.58442 | 2.898493 | 11.236520 | 0 |
| 52.88 | 97.03067 | 99.600374 | 25.026820 | 0 |
| 109.25 | 88.04000 | 0.000000 | 16.798451 | 0 |
| 105.60 | 104.09000 | 51.087736 | 6.731747 | 0 |
| 91.73 | 89.86937 | 160.752695 | 17.851526 | 0 |
| 91.97 | 95.35000 | 15.319511 | 14.058768 | 0 |
| 69.65 | 100.57019 | 76.930164 | 22.999587 | 0 |
| 111.50 | 58.99000 | 7.373436 | 8.195322 | 0 |
ggpubr::ggarrange(
Training_datos %>%
group_by(Var_respuesta) %>%
summarise(recuento = n()) %>%
ungroup() %>%
mutate(prop_var_respuesta = recuento / sum(recuento) ) %>%
ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
geom_bar(stat = "identity") +
geom_text(
aes(label = scales::percent(prop_var_respuesta)),
position = position_stack(vjust = 0.5), size = 4) +
scale_fill_manual( values = c("azure4", "azure3"))
,
Training_datos_adasyn %>%
group_by(Var_respuesta) %>%
summarise(recuento = n()) %>%
mutate(prop_var_respuesta = recuento / sum(recuento) ) %>%
ungroup() %>%
ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
geom_bar(stat = "identity") +
geom_text(
aes(label = scales::percent(prop_var_respuesta)),
position = position_stack(vjust = 0.5), size = 4) +
scale_fill_manual( values = c("azure4", "azure3")),
common.legend = T, legend = 'bottom'
,
labels = c('Normal', 'Sobre Muestreo\n Adasyn'), label.x = 0.1
) %>%
plot()Por el momento dejaré una receta lista para continuar con el taller.
# Creamos la Fórmula para la receta a partir
# de nombres de las variables y un poco de magia con paste
Receta_modelo_1_formula <- Training_datos %>%
dplyr::select(
dplyr::matches('Cog'), dplyr::matches('escala')) %>%
names() %>%
paste(., collapse = ' + ') %>%
paste('Var_respuesta ~' ,.) %>%
as.formula()
# 🔴 1.) Iniciamos la receta
Receta_modelo_1 <- recipe(
formula = Receta_modelo_1_formula,
data = Training_datos ,
strata = Var_respuesta) %>%
# 🔴 2.) Steps
# 🟡 2.1) Eliminamos las variables que contengan más de un 20% de datos perdidos
step_filter_missing(all_predictors(), threshold = 0.1) %>%
# 🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
step_impute_bag(all_numeric(),-all_outcomes()) %>%
# 2.3) Normalizamos los datos (restamos la media)
step_normalize(all_numeric(),-all_outcomes() ) %>%
# 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
step_scale(all_numeric(),-all_outcomes()) %>%
# 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>%
# 2.6) Sobre muestreamos la variable para equilibrar grupos
step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>%
# 3.) Preparamos la receta
prep()| Cognicion_01 | Escala_02 | Escala_03 | Var_respuesta |
|---|---|---|---|
| -0.6426710 | -0.0464912 | 1.5341196 | 0 |
| 0.7604689 | -0.7638606 | -0.3800329 | 0 |
| -2.9059967 | -0.8783295 | -1.4610980 | 0 |
| -1.1137735 | 0.7889501 | -0.7053657 | 0 |
| 0.7967336 | -0.2058789 | -1.3439820 | 0 |
| 0.6730268 | -1.4229671 | -0.9510280 | 0 |
| 0.2029411 | -0.0785598 | -0.0393877 | 0 |
| 0.2110752 | -0.5371131 | -0.1631135 | 0 |
| -0.5454002 | 0.5438529 | -0.4982340 | 0 |
| 0.8729912 | -1.2460174 | -1.4502475 | 0 |
Vale, dejamos atrás la parte del pre-procesamiento y empezamos la parte más machine leaninrg. Modelos. Modelos a la carrera.
La función del paquete parsnip dentro de tidymodels es propoconar una interfaz universal para diferentes métodos de modelado en R.
El ejemplo del principio del taller está basado en esta sección, parnsip es el encargado de que solo con cambiar unas pocas líneas puedas cambiar todo un modelo dentro de un flujo de aprendizaje automático.
A fecha de este taller, tidymodels cuenta con más de 155 modelos, incluyendo supervivencia y series de tiempo. Se pueden encontrar en: https://www.tidymodels.org/find/parsnip/
Supongamos que queremos empezar modelar, queremos usar una regresión logística, un clásico en este campo.Pero dentro de R hay varias librarías que usan la regresión logística, cada una con sus variaciones y reglas. Veamos que tenemos para esto:
## # A tibble: 8 × 2
## engine mode
## <chr> <chr>
## 1 glm classification
## 2 glmnet classification
## 3 LiblineaR classification
## 4 spark classification
## 5 keras classification
## 6 stan classification
## 7 brulee classification
## 8 h2o classification
logistic_reg(
engine = "glmnet",
#Hyper-parametros
penalty = NULL, mixture = NULL, mode = 'classification')## Logistic Regression Model Specification (classification)
##
## Computational engine: glmnet
logistic_reg(
engine = "glmnet",
#Hyper-parametros
penalty = 0, mixture = 0, mode = 'classification')## Logistic Regression Model Specification (classification)
##
## Main Arguments:
## penalty = 0
## mixture = 0
##
## Computational engine: glmnet
logistic_reg(
engine = "glmnet",
#Hyper-parametros
penalty = tune(), mixture = tune(), mode = 'classification')## Logistic Regression Model Specification (classification)
##
## Main Arguments:
## penalty = tune()
## mixture = tune()
##
## Computational engine: glmnet
Model_RandomForest <-
# 🔴 Especificamos el modelo que queremos en este caso un random forest
rand_forest(
# 🟡 Ajustamos los hiper-parametros
mtry = tune(), trees = tune(), min_n = tune()) %>%
# 🔴 Ponemos el motor, es decir la librería por la queremos que se ejecute el modelo
set_engine("ranger", importance = "impurity") %>%
# 🔴 Ajustamos el modo del modelo, es decri si queremos regresion o clasifiacion
set_mode("classification")El paquete de workflows permite combinar una receta y
una especificación de modelo en un único objeto. Esto hace que todo el
proceso de modelado, desde el pre-procesamiento de datos hasta el ajuste
del modelo, esté en un mismo punto
yardstickEl paquete yardstick, en concreto, se utiliza para las
métricas de evaluación de modelos. Además
# 🔴 Seleccionamos las métricas que queremos para nuestro modelo
Modelo_Metricas <- metric_set(accuracy, j_index, precision, sensitivity, specificity, roc_auc, f_meas,recall, mcc)También decir que es es posible crear tu propia métrica. Pero aún no he indagado en esta metodología. Se puede encontrar más en: https://yardstick.tidymodels.org/articles/custom-metrics.html
Los hiperparámetros son una pieza crucial en los modelos de aprendizaje automático. Estos son parámetros de configuración externos que no pueden aprenderse a partir de los datos, pero que influyen significativamente en el rendimiento del modelo. Un ejemplo sería el número de árboles de un random forest por ejemplo. A diferencia de los parámetros del modelo, que se aprenden durante el entrenamiento, los híperparámetros se establecen antes de que comience el entrenamiento de este.
La selección de los hiperparámetros adecuados es esencial para conseguir un modelo de aprendizaje automático generalizable y con buen rendimiento y es un proceso delicado. Los híper-parámetros desempeñan un papel crucial a la hora de controlar el equilibrio entre la sobreajuste y el infraajuste (overfitting y underfitting). El sobreajuste se produce cuando un modelo funciona bien con los datos de entrenamiento pero no consigue generalizarse a datos nuevos que no se han visto. Por otro lado, el infraajuste se produce cuando el modelo es demasiado simple y no puede captar los patrones subyacentes en los datos. Ajustar los hiperparámetros puede ayudar a encontrar el nivel adecuado de complejidad del modelo, haciendo un buen equilibrio entre la varianza del conjunto de entrenamiento y el de evaluación.
Los paquetes tune y dials son los que nos
permiten hacer esta selección y control. Aquí se crea todo un subflujo
de trabajo para el ajuste de híper-parámetros.
grid_regular(): hará todas las combinaciones entre los
rangos de la variable.grid_random(): probará hiper parámetros
aleatoriamente.grid_max_entropy(): propondrá una combinacion de
hiperpárametros que garanticen que se cubrirá todo el espectro, con fin
de evitar mínimos locales en la conversión de algoritmos.grid_latin_hypercube(): propondrá una matriz.grid_regular() hará las combinaciones acorde al rango
que le hayamos dado. Por defecto hace la combinación entre el mínimo, el
medio y el máximo de cada valor.
## # A tibble: 27 × 3
## mtry min_n trees
## <int> <int> <int>
## 1 2 2 2
## 2 6 2 2
## 3 10 2 2
## 4 2 6 2
## 5 6 6 2
## 6 10 6 2
## 7 2 10 2
## 8 6 10 2
## 9 10 10 2
## 10 2 2 6
## # ℹ 17 more rows
Podemos ampliar esto con el parámetro levels. Ahora será con el mínimo, el quantile 33, el cuantil 66 y el máximo.
grid_regular(
mtry(range = c(2, 10)),
min_n(range = c(2, 10)),
trees(range = c(2, 10)),
levels = 4 )## # A tibble: 64 × 3
## mtry min_n trees
## <int> <int> <int>
## 1 2 2 2
## 2 4 2 2
## 3 7 2 2
## 4 10 2 2
## 5 2 4 2
## 6 4 4 2
## 7 7 4 2
## 8 10 4 2
## 9 2 7 2
## 10 4 7 2
## # ℹ 54 more rows
grid_max_entropy(
mtry(range = c(1, 4)),
min_n(range = c(10, 30)),
trees(range = c(1, 1000)),
size = 10
)## # A tibble: 10 × 3
## mtry min_n trees
## <int> <int> <int>
## 1 1 25 20
## 2 2 11 610
## 3 4 12 2
## 4 2 20 936
## 5 4 12 733
## 6 4 30 868
## 7 1 20 676
## 8 2 19 307
## 9 2 29 661
## 10 4 26 190
yardstick con dial
y tuneWF_hiper_parametros <-
# 🔴 Este proceso se hace con la función tune grid 🔴
tune_grid(
# Ponemos la receta nuestro modelo, que incuye el preprocesamiento y la fórmula
object = WF_Receta_modelo_1_Random_forest,
# 🟡 Ponemos los Folds del conjunto de entrenamiento,
# 🟡 Dentro del entrenamiento hará un sub entrenamiento y una sub validación
resamples = Folds_training,
#Ponemos un grid para que pruebe con diferentes híper-parámetros, ya con un rango pre-definido
grid = grid_max_entropy(
mtry(range = c(1, 4)),
min_n(range = c(10, 30)),
trees(range = c(1, 1000)),
size = 10),
# Ponemos las métricas que hemos definido anteriormente
metrics = Modelo_Metricas,
# Con esta opción podemos guardar las predicciones
control = control_grid( save_pred = T)
)WF_hiperparametros_coleccion <-
WF_hiper_parametros %>%
# 🔴 Esta es la funcion que retorna las métricas, importante tenerla en cuenta
collect_metrics()
WF_hiperparametros_coleccion %>%
head(.,20) %>%
knitr::kable()| mtry | trees | min_n | .metric | .estimator | mean | n | std_err | .config |
|---|---|---|---|---|---|---|---|---|
| 2 | 886 | 10 | accuracy | binary | 0.6610256 | 5 | 0.0496335 | Preprocessor1_Model01 |
| 2 | 886 | 10 | f_meas | binary | 0.7673926 | 5 | 0.0365573 | Preprocessor1_Model01 |
| 2 | 886 | 10 | j_index | binary | 0.1832661 | 5 | 0.1253188 | Preprocessor1_Model01 |
| 2 | 886 | 10 | mcc | binary | 0.1626371 | 5 | 0.1086564 | Preprocessor1_Model01 |
| 2 | 886 | 10 | precision | binary | 0.8400397 | 5 | 0.0310967 | Preprocessor1_Model01 |
| 2 | 886 | 10 | recall | binary | 0.7082661 | 5 | 0.0439273 | Preprocessor1_Model01 |
| 2 | 886 | 10 | roc_auc | binary | 0.7213836 | 5 | 0.0469202 | Preprocessor1_Model01 |
| 2 | 886 | 10 | sensitivity | binary | 0.7082661 | 5 | 0.0439273 | Preprocessor1_Model01 |
| 2 | 886 | 10 | specificity | binary | 0.4750000 | 5 | 0.0918559 | Preprocessor1_Model01 |
| 4 | 123 | 28 | accuracy | binary | 0.6760256 | 5 | 0.0501729 | Preprocessor1_Model02 |
| 4 | 123 | 28 | f_meas | binary | 0.7747366 | 5 | 0.0350463 | Preprocessor1_Model02 |
| 4 | 123 | 28 | j_index | binary | 0.2768145 | 5 | 0.1524681 | Preprocessor1_Model02 |
| 4 | 123 | 28 | mcc | binary | 0.2346557 | 5 | 0.1298836 | Preprocessor1_Model02 |
| 4 | 123 | 28 | precision | binary | 0.8684292 | 5 | 0.0392609 | Preprocessor1_Model02 |
| 4 | 123 | 28 | recall | binary | 0.7018145 | 5 | 0.0393757 | Preprocessor1_Model02 |
| 4 | 123 | 28 | roc_auc | binary | 0.7346900 | 5 | 0.0462352 | Preprocessor1_Model02 |
| 4 | 123 | 28 | sensitivity | binary | 0.7018145 | 5 | 0.0393757 | Preprocessor1_Model02 |
| 4 | 123 | 28 | specificity | binary | 0.5750000 | 5 | 0.1286954 | Preprocessor1_Model02 |
| 2 | 81 | 16 | accuracy | binary | 0.6762821 | 5 | 0.0479575 | Preprocessor1_Model03 |
| 2 | 81 | 16 | f_meas | binary | 0.7791279 | 5 | 0.0342760 | Preprocessor1_Model03 |
WF_hiper_parametros %>%
collect_metrics() %>%
pivot_longer(cols = c('mtry','trees', 'min_n'), names_to = 'tipo_parametro', values_to = 'valor_parametro') %>%
ggplot(aes(valor_parametro,mean, color = mean)) +
geom_point() +
facet_grid(.metric~tipo_parametro, scales = 'free' ) +
scale_color_viridis_b() +
theme_bw()WF_hiperparametros_mejor <- WF_hiper_parametros %>% show_best("roc_auc",1)
WF_hiperparametros_mejor %>%
knitr::kable()| mtry | trees | min_n | .metric | .estimator | mean | n | std_err | .config |
|---|---|---|---|---|---|---|---|---|
| 2 | 240 | 29 | roc_auc | binary | 0.7387727 | 5 | 0.0384544 | Preprocessor1_Model06 |
fit y predictYa hemos hecho todos los procesos previos. Todo lo que a un modelo de machine learning le puede hacer falta. Finalmente, una vez que lo tenemos todo hilado, hay que proceder a entrenar el modelo y dejarlo listo.
Modelo_fitted <- WF_Receta_modelo_1_Random_forest %>%
# 🔴 Finalizamos el flujo de trabajo con la mejor combinación de Hiper paraemtros encontrada
finalize_workflow(WF_hiperparametros_mejor) %>%
# 🔴 finalmente ajustamos con fit
fit(data = Training_datos)Con fit ya por fin tenemos nuestro modelo, ya podemos
ponernos a hacer lo que hace el aprendizaje automático: predecir. vamos
a crear una observación nueva y a ver como la trabaja.
Nuevo_caso <- tibble(
'Cognicion_01' = 123,
'Cognicion_02' = 102,
'Cognicion_03' = 89,
'Cognicion_04' = 100,
'Cognicion_05' = 179,
'Escala_01' = 12,
'Escala_02' = 22,
'Escala_03' = 20
)tibble(
# 🔴 La función predict de la librería stats tiene sinergias con tidymodels, neecsitra de estos argumentos:
# - 🟡 El modelo ya ajustado (model_fitted)
# - 🟡 Un tibble/dataframe con la nueva prediccion (debe contener las mismas variables)
# - 🟡 Un tipo de prediccion en los casos de clasificación, si queremos las probabiliadades de cada categoria o la categoría predicha directamente
predict(Modelo_fitted, Nuevo_caso, type = "prob" ),
predict(Modelo_fitted, Nuevo_caso, type = "class")) %>%
knitr::kable()| .pred_0 | .pred_1 | .pred_class |
|---|---|---|
| 0.9541702 | 0.0458298 | 0 |
tibble(predict(Modelo_fitted, Testing_datos, type = "prob" ), predict(Modelo_fitted, Testing_datos, type = "class")) %>%
head(.,10) %>%
knitr::kable()| .pred_0 | .pred_1 | .pred_class |
|---|---|---|
| 0.7868158 | 0.2131842 | 0 |
| 0.5619997 | 0.4380003 | 0 |
| 0.8653084 | 0.1346916 | 0 |
| 0.4211667 | 0.5788333 | 1 |
| 0.8627070 | 0.1372930 | 0 |
| 0.9247541 | 0.0752459 | 0 |
| 0.0893169 | 0.9106831 | 1 |
| 0.9033712 | 0.0966288 | 0 |
| 0.9443273 | 0.0556727 | 0 |
| 0.7467458 | 0.2532542 | 0 |
last_fitCrear las predicciones y todo es un sistema engorroso. Además de que se pierde mucha informaciíon por el camino, ya que todo sigue estando en objetos dispersos. Para ello los desarrolladores crearon last_fit.
La función last_fit permite tener todo unificado en un solo tibble, al tener que ajustarse con un objeto split, esta recibe directamente los datos de entrenamiento y los datos de evaluación (recibe datos de training y testing) y entonces crea sus métricas y sus predicciones.
Con last_fit() contamos con una serie de funciones para
extraer todo lo que necesitemos de ellas:
extract_fit_engine() : permite extraer el modelo crudo
como si no hubiera sido ajustado con tidymodels, de modo que si es un
modelo explicativo puede ser tratado con normalidadcollect_metrics() : extraer un tibble con las métricas
del modelo (puestas con metric set)collect_predictions() : extrae un tibble con las
predicciones de los datos de evaluación. (Por defecto el umbral de
rendimiento en esta función es 0.5, y estoy buscando alguna forma de
cambiar esto, pero se puede manipular el mismo tibble como cualquier
otro, de modo que no es un drama.)extract_workflow() : permite extraer todo el flujo del
modelo, de incluyendo la receta de modo que podemos hacer
predicciones:Modelo_1_final <- WF_Receta_modelo_1_Random_forest %>%
# 🔴 Finalizamos el flujo de trabajo con la mejor combinación de Hiper paraemtros encontrada
finalize_workflow(WF_hiperparametros_mejor) %>%
# 🔴 Ajustamos con last_fit, poniendo el split de datos inicial y las metricas para las predicciones
last_fit(Split_datos_strat, metrics = Modelo_Metricas)
Modelo_1_final## # Resampling results
## # Manual resampling
## # A tibble: 1 × 6
## splits id .metrics .notes .predictions .workflow
## <list> <chr> <list> <list> <list> <list>
## 1 <split [198/87]> train/test split <tibble> <tibble> <tibble> <workflow>
Modelo_1_final %>%
# 🔴 dentro del flujo de trabajo están la receta y el modelo ya ajustado
extract_workflow() %>%
# de modo que lo podemos usar para hacer predicciones
predict(Testing_datos, type = "prob" )## # A tibble: 87 × 2
## .pred_0 .pred_1
## <dbl> <dbl>
## 1 0.787 0.213
## 2 0.562 0.438
## 3 0.865 0.135
## 4 0.421 0.579
## 5 0.863 0.137
## 6 0.925 0.0752
## 7 0.0893 0.911
## 8 0.903 0.0966
## 9 0.944 0.0557
## 10 0.747 0.253
## # ℹ 77 more rows
A partir de aquí podemos manipular el objeto fit de múltiples maneras múltiples.
Modelo_1_final %>%
# 🔴 Con esta funcion podemos recopilar todas las metricas de un modelo
collect_metrics() %>%
select(-.estimator, -.config) %>%
mutate(.estimate = round(.estimate,2)) %>%
knitr::kable()| .metric | .estimate |
|---|---|
| accuracy | 0.71 |
| j_index | 0.39 |
| precision | 0.89 |
| sensitivity | 0.72 |
| specificity | 0.67 |
| f_meas | 0.80 |
| recall | 0.72 |
| mcc | 0.33 |
| roc_auc | 0.74 |
options(yardstick.event_first = FALSE)
Modelo_1_final %>% collect_predictions() %>%
# 🔴 Esta funcion permite calcular rapidamente la curva roc de un modelo ya ajustado
roc_curve(truth = Var_respuesta, .pred_0 ) %>%
ggplot(aes( x = 1 - specificity ,y = sensitivity ) ) +
geom_path() +
geom_abline(intercept = 0 , slope = 1, linetype = 3) +
labs(title = 'Cruva roc del modelo_1_final') +
theme_bw()Advertencia: en la función de curva roc, ved que
puesto predicción del 0 en vez del 1. Esto es un error interno en la
clasificación dentro del del paquete yardstick. Hay que
tener cuidado sobre si estas funciones utilizan el primer / segundo / x
nivel de su factor como el “evento”. Por defecto, yardstick
elige el primer nivel de verdad como “evento” al calcular la curva roc.
Podéis encontrar más de este suceso en el github de los desarrolladores:
https://github.com/tidymodels/yardstick/issues/94
El concepto de rendimiento del umbral (threshold performance) es crucial a la hora de evaluar y ajustar modelos de clasificación binaria.
El umbral es el valor de probabilidad por encima del cual un resultado predicho se considera de clase positiva. Por defecto, muchos modelos utilizan un umbral de 0,5. Sin embargo, en casi todos los casos, el umbral óptimo no estará otro punto. Este “óptimo lo tendremos que definir que en función de los objetivos y requisitos específicos del problema.
Ajustar el umbral permite elegir entre precisión y la recall. Bajar el umbral aumenta la sensibilidad (recall) pero disminuye la precisión, y viceversa. Básicamente, el umbral óptimo depende de la importancia relativa de los falsos positivos y los falsos negativos en una aplicación concreta.
En tidymodels está la librería probably con la función
threshold_perf(), la cual nos permite, dadas las
predicciones del modelo ver como se calcula este umbral.
pacman::p_load(probably)
Modelo_1_final_Threshold_performance <-
# Llamamos al modelo ya finallizado con "last_fit"
Modelo_1_final %>%
# Coleccionamos sus predicciones
collect_predictions() %>%
# 🔴 Ejecutamos el redimeinto de umbral con la siguiente función
threshold_perf(
# Le decimos que variable es la verdadera respuesta
truth = Var_respuesta,
# le decimos que variable es la predicción (aviso: debido a un bugg en la version )
.pred_1,
# ajustamos un rango en el que se evaluará el umbral en cada punto
thresholds = seq(0.5,1,0.01) )
# displayment
Modelo_1_final_Threshold_performance %>%
head(., 10) %>%
mutate(.estimate = round(.estimate, 2)) %>%
knitr::kable()| .threshold | .metric | .estimator | .estimate |
|---|---|---|---|
| 0.50 | sensitivity | binary | 0.28 |
| 0.51 | sensitivity | binary | 0.28 |
| 0.52 | sensitivity | binary | 0.28 |
| 0.53 | sensitivity | binary | 0.28 |
| 0.54 | sensitivity | binary | 0.28 |
| 0.55 | sensitivity | binary | 0.28 |
| 0.56 | sensitivity | binary | 0.26 |
| 0.57 | sensitivity | binary | 0.25 |
| 0.58 | sensitivity | binary | 0.23 |
| 0.59 | sensitivity | binary | 0.23 |
Modelo_1_final_Threshold_performance %>%
ggplot(aes(x = .threshold, y = .estimate, color = .metric)) +
geom_line(size = 1.1) +
scale_color_manual(values = c('purple','blue','pink')) +
geom_vline(xintercept =
Modelo_1_final_Threshold_performance %>%
filter(.metric == 'j_index') %>%
arrange(desc(.estimate)) %>%
slice(1) %>%
pull(.threshold), linetype = 'dashed', size = 1.5)workflow_set()Bien, ya hemos aprendido a crear todo un flujo entero de machine learning con tidymodels, desde la base hasta su final. Muy útil todo. pero calibrar modelos de 1 en 1 es un auténtico petardazo. No me acaba de ser útil este taller. ¿Cómo me va a ser útil si tengo que probar unas 20 fórmulas con datos?
Lo entiendo perfectamente, a mi también me pasó. Los médicos me preguntaron que era mejor para tener en cuenta, si las escalas, las cogniciones, la bioquímica o la genética a la hora de diagnosticar al paciente.
Comprobémoslo…
Receta_modelo_Escalas_formula <- Training_datos %>%
dplyr::select(dplyr::matches('Esc')) %>%
names() %>%
paste(., collapse = ' + ') %>%
paste('Var_respuesta ~' ,.) %>%
as.formula()
# 🔴 1.) Iniciamos la receta
Receta_modelo_Escalas <- recipe(
formula = Receta_modelo_Escalas_formula,
data = Training_datos ,
strata = Var_respuesta) %>%
# 🔴 2.) Steps
# 🟡 2.1) Eliminamos las variables que contengan más de un 25% de datos perdidos
step_filter_missing(all_predictors(), threshold = 0.25) %>%
# 🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
step_impute_bag(all_numeric(),-all_outcomes()) %>%
# 2.3) Normalizamos los datos (restamos la media)
step_normalize(all_numeric(),-all_outcomes() ) %>%
# 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
step_scale(all_numeric(),-all_outcomes()) %>%
# 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>%
# 2.6) Sobre muestreamos la variable para equilibrar grupos
step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>%
# 3.) Preparamos la receta
prep()
Receta_modelo_Cognicion_formula <- Training_datos %>%
dplyr::select(
dplyr::matches('Cog')) %>%
names() %>%
paste(., collapse = ' + ') %>%
paste('Var_respuesta ~' ,.) %>%
as.formula()
# 🔴 1.) Iniciamos la receta
Receta_modelo_Cognicion <- recipe(
formula = Receta_modelo_Cognicion_formula,
data = Training_datos ,
strata = Var_respuesta) %>%
# 🔴 2.) Steps
# 🟡 2.1) Eliminamos las variables que contengan más de un 25% de datos perdidos
step_filter_missing(all_predictors(), threshold = 0.25) %>%
# 🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
step_impute_mean(all_numeric(),-all_outcomes()) %>%
# 2.3) Normalizamos los datos (restamos la media)
step_normalize(all_numeric(),-all_outcomes() ) %>%
# 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
step_scale(all_numeric(),-all_outcomes()) %>%
# 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>%
# 2.6) Sobre muestreamos la variable para equilibrar grupos
step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>%
# 3.) Preparamos la receta
prep()
bake(Receta_modelo_Cognicion, new_data = NULL)## # A tibble: 316 × 4
## Cognicion_01 Cognicion_03 Cognicion_04 Var_respuesta
## <dbl> <dbl> <dbl> <fct>
## 1 -0.643 -2.70e-15 2.51e-15 0
## 2 0.760 3.13e- 1 1.13e+ 0 0
## 3 -2.91 -2.70e-15 2.51e-15 0
## 4 -1.11 -2.70e-15 2.51e-15 0
## 5 0.797 1.64e+ 0 2.59e- 1 0
## 6 0.673 5.70e- 1 3.93e- 1 0
## 7 0.203 -2.70e-15 2.51e-15 0
## 8 0.211 2.03e- 1 8.86e- 1 0
## 9 -0.545 -2.70e-15 2.51e-15 0
## 10 0.873 1.50e+ 0 1.09e+ 0 0
## # ℹ 306 more rows
Receta_modelo_Bioquimicas_formula <- Training_datos %>%
dplyr::select(
dplyr::matches('Bio')) %>%
names() %>%
paste(., collapse = ' + ') %>%
paste('Var_respuesta ~' ,.) %>%
as.formula()
# 🔴 1.) Iniciamos la receta
Receta_modelo_Bioquimica <- recipe(
formula = Receta_modelo_Bioquimicas_formula,
data = Training_datos ,
strata = Var_respuesta) %>%
# 🔴 2.) Steps
# 🟡 2.1) Eliminamos las variables que contengan más de un 25% de datos perdidos
step_filter_missing(all_predictors(), threshold = 0.25) %>%
# 🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
step_impute_mean(all_numeric()) %>%
# 2.3) Normalizamos los datos (restamos la media)
step_normalize(all_numeric(),-all_outcomes() ) %>%
# 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
step_scale(all_numeric(),-all_outcomes()) %>%
# 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>%
# 2.6) Sobre muestreamos la variable para equilibrar grupos
step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>%
# 3.) Preparamos la receta
prep()
Receta_modelo_Popurri <- Training_datos %>%
dplyr::select(
Cognicion_01, Cognicion_02, Escala_01,Escala_02, Bioquimica_01,Bioquimica_02, Genetica_01,Genetica_02) %>%
names() %>%
paste(., collapse = ' + ') %>%
paste('Var_respuesta ~' ,.) %>%
as.formula()
# 🔴 1.) Iniciamos la receta
Receta_modelo_Popurri <- recipe(
formula = Receta_modelo_Popurri,
data = Training_datos ,
strata = Var_respuesta) %>%
# 🔴 2.) Steps
# 🟡 2.1) Eliminamos las variables que contengan más de un 5% de datos perdidos
step_filter_missing(all_predictors(), threshold = 0.25) %>%
# 🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
step_impute_bag(all_numeric(),-all_outcomes()) %>%
# 2.3) Normalizamos los datos (restamos la media)
step_normalize(all_numeric(),-all_outcomes() ) %>%
# 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
step_scale(all_numeric(),-all_outcomes()) %>%
# 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>%
# 2.6) Sobre muestreamos la variable para equilibrar grupos
step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>%
# 3.) Preparamos la receta
prep()workflow_set() para unificarlas todasAhora, esto son muchas recetas y mucho texto. Será difícil combinarlo
todo, ¿no?…¡Qué va!. workflow_set() ya está pensado para
esto. Funciona como el crossing() de la librería purrr o
como el expand.grid() de la librería base,
crear workflows con todas las combinaciones.
Receta_list <- list(
'Modelo_Escalas' = Receta_modelo_Escalas,
'Modelo_Cognicion' = Receta_modelo_Cognicion,
'Modelo_Bioquimica' = Receta_modelo_Bioquimica,
'Modelo_Popurri' = Receta_modelo_Popurri
)
Model_list <- list(
RandomForest = Model_RandomForest
# ,XGBoost = Model_XGBoost
# ,Bagged_Trees = Model_Bagged_Tree
)
Wokflows_set <- workflow_set(
preproc = Receta_list,
models = Model_list,
# 🔴 La opción Cross es para que haga todos los cruces de modelos con recetas
cross = T)
Wokflows_set## # A workflow set/tibble: 4 × 4
## wflow_id info option result
## <chr> <list> <list> <list>
## 1 Modelo_Escalas_RandomForest <tibble [1 × 4]> <opts[0]> <list [0]>
## 2 Modelo_Cognicion_RandomForest <tibble [1 × 4]> <opts[0]> <list [0]>
## 3 Modelo_Bioquimica_RandomForest <tibble [1 × 4]> <opts[0]> <list [0]>
## 4 Modelo_Popurri_RandomForest <tibble [1 × 4]> <opts[0]> <list [0]>
Una vez hecho el workflow_set, se pone a trabajar con
workflow_map.
Wokflows_set_map <-
Wokflows_set %>%
# 🔴 workflow maps es una funcion que permite ajustar múltiples flujos de trabajo
workflow_map(
resamples = Folds_training,
fn = "tune_grid",
grid = grid_max_entropy(
mtry(range = c(1, 4)),
min_n(range = c(10, 30)),
trees(range = c(1, 1000)),
size = 10),
# verbose = TRUE,
metrics = Modelo_Metricas,
control = control_grid( save_pred = T),
# 🔴 para garantizare replicabildiad, podemos fijar la semilla en el proceso
seed = 2465)
Wokflows_set_map## # A workflow set/tibble: 4 × 4
## wflow_id info option result
## <chr> <list> <list> <list>
## 1 Modelo_Escalas_RandomForest <tibble [1 × 4]> <opts[4]> <tune[+]>
## 2 Modelo_Cognicion_RandomForest <tibble [1 × 4]> <opts[4]> <tune[+]>
## 3 Modelo_Bioquimica_RandomForest <tibble [1 × 4]> <opts[4]> <tune[+]>
## 4 Modelo_Popurri_RandomForest <tibble [1 × 4]> <opts[4]> <tune[+]>
AVISO, workflow_map puede
paralelizarse, pero aún estoy descubriendo como integrarlo de forma
efectiva.
Finalmente, podemos calibrar cada modelo individualmente y tener todo finiquitado. Ese proceso sería más eficiente con una función, desde luego. Pero se me acaba el tiempo para dedicar al taller :(
# Modelo Escalas ----
Modelo_Escalas_RandomForest_hyperparametros <- Wokflows_set_map %>%
# 🔴 Extraemos el flujo correspondiente al modelo que queremos
extract_workflow_set_result('Modelo_Escalas_RandomForest') %>%
# 🔴 podemos ver una tabla con los mejores resultados, acorde a una metrica
show_best(metric = 'sensitivity', n = 50) %>%
filter(row_number() == 1) %>%
select(mtry,trees,min_n,.config )
Modelo_Escalas_RandomForest_final <- Wokflows_set_map %>%
# 🔴 Extraemos el flujo correspondiente al modelo que queremos
extract_workflow('Modelo_Escalas_RandomForest') %>%
finalize_workflow(Modelo_Escalas_RandomForest_hyperparametros) %>%
last_fit(Split_datos_strat, metrics = Modelo_Metricas )
# Modelo_Cognicion ----
Modelo_Cognicion_RandomForest_hyperparametros <- Wokflows_set_map %>%
extract_workflow_set_result('Modelo_Cognicion_RandomForest') %>%
show_best(metric = 'sensitivity', n = 50) %>%
filter(row_number() == 1) %>%
select(mtry, trees, min_n, .config )
Modelo_Cognicion_RandomForest_final <- Wokflows_set_map %>%
extract_workflow('Modelo_Cognicion_RandomForest') %>%
finalize_workflow(Modelo_Cognicion_RandomForest_hyperparametros) %>%
last_fit(Split_datos_strat, metrics = Modelo_Metricas )
# Modelo_Bioquimico ----
Modelo_Bioquimica_RandomForest_hyperparametros <- Wokflows_set_map %>%
extract_workflow_set_result('Modelo_Bioquimica_RandomForest') %>%
show_best(metric = 'sensitivity', n = 50) %>%
filter(row_number() == 1) %>%
select(mtry, trees, min_n, .config )
Modelo_Bioquimica_RandomForest_final <- Wokflows_set_map %>%
extract_workflow('Modelo_Bioquimica_RandomForest') %>%
finalize_workflow(Modelo_Bioquimica_RandomForest_hyperparametros) %>%
last_fit(Split_datos_strat, metrics = Modelo_Metricas )
# Modelo Popurri ----
Modelo_Popurri_RandomForest_hyperparametros <- Wokflows_set_map %>%
extract_workflow_set_result('Modelo_Popurri_RandomForest') %>%
show_best(metric = 'sensitivity', n = 50) %>%
filter(row_number() == 1) %>%
select(mtry, trees, min_n, .config )
Modelo_Popurri_RandomForest_final <- Wokflows_set_map %>%
extract_workflow('Modelo_Popurri_RandomForest') %>%
finalize_workflow(Modelo_Popurri_RandomForest_hyperparametros) %>%
last_fit(Split_datos_strat, metrics = Modelo_Metricas )Una vez acabado todo, podemos ver todo lo necesario usando unos pocos
verbos como collect() y las lógicas del tidyverse.
map2(
list(
Modelo_Escalas_RandomForest_final,
Modelo_Cognicion_RandomForest_final,
Modelo_Bioquimica_RandomForest_final,
Modelo_Popurri_RandomForest_final),
list('Modelo_Escala','Modelo_Cognicion','Modelo_Bioquimica','Modelo_Popurri'),
~ ..1 %>%
collect_metrics() %>%
mutate('Modelo'= ..2) %>%
select(.metric,.estimate,Modelo )) %>%
bind_rows() %>%
pivot_wider(names_from = Modelo, values_from = .estimate) %>%
knitr::kable()| .metric | Modelo_Escala | Modelo_Cognicion | Modelo_Bioquimica | Modelo_Popurri |
|---|---|---|---|---|
| accuracy | 0.6436782 | 0.6091954 | 0.7471264 | 0.6666667 |
| j_index | 0.2222222 | -0.0265700 | 0.0241546 | 0.2101449 |
| precision | 0.8518519 | 0.7868852 | 0.7974684 | 0.8448276 |
| sensitivity | 0.6666667 | 0.6956522 | 0.9130435 | 0.7101449 |
| specificity | 0.5555556 | 0.2777778 | 0.1111111 | 0.5000000 |
| f_meas | 0.7479675 | 0.7384615 | 0.8513514 | 0.7716535 |
| recall | 0.6666667 | 0.6956522 | 0.9130435 | 0.7101449 |
| mcc | 0.1855216 | -0.0235126 | 0.0338612 | 0.1805788 |
| roc_auc | 0.7504026 | 0.5024155 | 0.5048309 | 0.6231884 |
ROC_CURVE_training_todos_los_posibles_modelos <- Wokflows_set_map %>%
collect_predictions() %>%
group_split(wflow_id) %>%
set_names(list('Escala','Cognicion','Bioquimica','Popurri')) %>%
map(~ .x %>% roc_curve(truth = Var_respuesta , .pred_0)) %>%
map2(., list('Escala','Cognicion','Bioquimica','Popurri'),
~ .x %>% mutate(model = .y)) %>%
bind_rows()
Plot_ROC_CURVE_training_todos_los_posibles_modelos <- ROC_CURVE_training_todos_los_posibles_modelos %>%
ggplot(aes( x = 1 - specificity,y = sensitivity, color = model ) ) +
geom_path() +
geom_abline(intercept = 0, slope = 1, linetype = 3) +
theme_bw()
ROC_CURVE_training_modelos_finales <- list(
list('Esc' ,Modelo_Escalas_RandomForest_hyperparametros$.config),
list('Cog' ,Modelo_Cognicion_RandomForest_hyperparametros$.config),
list('Bio' ,Modelo_Bioquimica_RandomForest_hyperparametros$.config),
list('Pop' ,Modelo_Popurri_RandomForest_hyperparametros$.config)) %>%
map(~ Wokflows_set_map %>%
collect_predictions() %>%
filter(
str_detect(wflow_id, .x[[1]]) &
.config== .x[[2]]
)) %>%
map(~ .x %>%
mutate(Model= str_extract(wflow_id, '^.{2}') ) %>%
select(Model,Var_respuesta, .pred_0, .pred_1, .pred_class)
) %>%
set_names(list('Escala','Cognicion','Bioquimica','Popurri') ) %>%
map(~ .x %>% roc_curve(truth = Var_respuesta, .pred_0)) %>%
map2(., list('Escala','Cognicion','Bioquimica','Popurri'),
~ .x %>% mutate(model= .y)) %>%
bind_rows()
Plot_ROC_CURVE_training_modelos_finales <- ROC_CURVE_training_modelos_finales %>%
ggplot(aes( x = 1 - specificity ,y = sensitivity, color = model ) ) +
geom_path() +
geom_abline(intercept = 0, slope = 1, linetype = 3) +
theme_bw()
ROC_CURVE_testing_modelos_finales <- list(
Modelo_Escalas_RandomForest_final,
Modelo_Cognicion_RandomForest_final,
Modelo_Bioquimica_RandomForest_final,
Modelo_Popurri_RandomForest_final
) %>%
map(~.x %>%
collect_predictions() %>%
roc_curve(Var_respuesta, .pred_0)) %>%
map2(., list('Escala','Cognicion','Bioquimica','Popurri'),
~ .x %>% mutate(model = .y)) %>%
bind_rows()
Plot_ROC_CURVE_testing_modelos_finales <- ROC_CURVE_testing_modelos_finales %>%
ggplot(aes(x = 1 - specificity ,y = sensitivity, color = model ) ) +
geom_path() +
geom_abline(intercept = 0, slope = 1, linetype = 3) +
theme_bw()ggpubr::ggarrange( Plot_ROC_CURVE_training_todos_los_posibles_modelos,
Plot_ROC_CURVE_training_modelos_finales,
Plot_ROC_CURVE_testing_modelos_finales ,
nrow = 1 ,
common.legend = T,
legend = "bottom",
labels = c('Todos_training','Training','Test')) Los valores de Shapley (shap values), que deben su nombre al premio Nobel Lloyd Shapley, son un concepto de la teoría de juegos cooperativos que ha acabado encontrado aplicaciones en el aprendizaje automático.
En teoría de juegos, los valores de Shapley son la distribución equitativa de la contribución de cada jugador en un juego cooperativo entre todos los jugadores. En el contexto del aprendizaje automático, las características o variables se consideran jugadores, y el juego consiste en cuánto contribuye cada característica a la predicción de un modelo.
Es decir, los valores de Shapley son sencillamente las marginales de cada variable dentro del modelo. Compilarlos es computacionalmente extenso, ya que hay que hacer cada combinación de valores para una fila y evaluar de cuanto es su marginal. El valor final llamado Shap value es promedio de todas las aportaciones de una variable marginal.
Aparte de su extenso coste computacional, todo shap value asume que está utilizando un modelo de caja negra. de manera que si el modelo no está bien especificado en tema de interacciones, etc. es posible que estén dando indicaciones erróneas.
Actualmente se pueden compilar en tidymodels utilizando la librería
agua, al ser una extensión de tidymodels del entorno de
h2o. El problema es que esto hace que el motor de los
modelos deba ser a la fuerza h2o.
Una propuesta sería buscar un modelo bien calibrado con tidymodels y pasarlo por h2o para hacer un diagnostico más eficaz.
h2o_start()
h2o::h2o.no_progress()
modelo_h20 <- rand_forest() %>%
set_engine("h2o", max_runtime_secs = 20) %>%
set_mode('classification') %>%
fit(Var_respuesta ~ ., data = Receta_modelo_Bioquimica %>% bake(new_data = NULL) )
entrada_h20 <- h2o::as.h2o(
Receta_modelo_Bioquimica %>% bake(new_data = NULL) %>%
dplyr::rename(.outcome = Var_respuesta))
modelo_h20 %>%
extract_fit_engine() %>%
h2o::h2o.shap_summary_plot(entrada_h20)Para cada predicción, los valores de Shapley representan la contribución de cada característica a la diferencia entre la predicción del modelo y la predicción media. En otras palabras, indican cuánto añade o resta cada característica a la predicción media.
Si una predicción media es de 0.5, y el shap value es de 0.1, indicará que esa variable sola podría llegar a tener un impacto tal que haga que la predicción llegue a 0.6.